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Abstract. This paper deals with the formulation and numerical implementation of a fully cou- 
pled continuum model for deformation-diffusion in linearized elastic solids. The mathematical 
model takes into account the effect of the deformation on the diffusion process, and the affect of 
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, the transport of an inert chemical species on the deformation of the solid. We then present a 

(— i . robust computational framework for solving the proposed mathematical model, which consists of 

coupled non-linear partial differential equations. It should be noted that many popular numerical 
formulations may produce unphysical negative values for the concentration, particularly, when the 
diffusion process is anisotropic. The violation of the non-negative constraint by these numerical 
formulations is not mere numerical noise. In the proposed computational framework we employ a 
' novel numerical formulation that will ensure that the concentration of the diffusant be always non- 

negative, which is one of the main contributions of this paper. Representative numerical examples 
are presented to show the robustness, convergence, and performance of the proposed computational 
framework. Another contribution of this paper is to systematically study the affect of transport 
of the diffusant on the deformation of the solid and vice-versa, and their implication in modeling 
, degradation/healing of materials. We show that the coupled response is both qualitatively and 

^SJ , quantitatively different from the uncoupled response. 
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1. INTRODUCTION 

In this paper we present a continuum mathematical model and a computational framework for 
degradation/healing in elastic solids due to the presence of a solute. We shall neglect chemical 
reactions as well as thermal effects. We shall also assume that the strains in the solid are small. 
The deformation is coupled with the diffusion process, and the diffusion process is in turn coupled 
with the deformation of the solid. The chosen problem belongs to a broader class of problems, 
namely, coupled deformation-diffusion problems. 

Coupled deformation-diffusion problems arise in many civil engineering, material science, and 
polymer science applications. An important example is degradation/healing of materials and struc- 
tures. Many man-made and natural materials degrade/heal due to environmental conditions, and 
structural components and superstructures are constantly exposed to adverse conditions. The fate 
of the transport of a diffusant will in turn depend on the deformation of the solid. Some other 
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examples of coupled deformation-diffusion problems are moisture damage in cementitious materials 
[501 EI] and asphalt [TTJ E2] , hydrogen embrittlement [EEJ EZ] , curing of ceramics [35] , and swelling 
of polymers and composites [581 155j. 

The governing equations of a deformation-diffusion problem are coupled non-linear partial dif- 
ferential equations. That is, the mechanical variables (strains and stresses) and diffusive variables 
(concentration and mass flux) are coupled through constitutive relations. It is (in general) not pos- 
sible to obtain analytical solutions for these kinds of problems, and one has to resort to numerical 
techniques to solve practical problems. Despite the importance of coupled deformation-diffusion 
problems, there is no robust and reliable computational framework for solving such coupled prob- 
lems. In particular, the existing numerical studies have one or more of the following limitations: 

• considered academic and unrealistic problems like infinite slabs and infinite cylinders, 

• did not consider anisotropic diffusion and/or assumed the medium to be homogeneous, or 

• did not consider the fact that conventional numerical formulations and finite element pack- 
ages produce unphysical negative values for the concentration in solving diffusion- type equa- 
tions. 

1.1. Maximum principles and non- negative solutions for diffusion-type equations. Pre- 
dictive numerical simulations require accurate and reliable discretization methods. The resulting 
discrete systems must inherit or mimic fundamental properties of continuous systems. Maximum 
principles form an important set of properties for diffusion-type equations as these maximum princi- 
ples have mathematical implications and physical consequences. In the study on partial differential 
equations, maximum principles are often used in existence theorems, and in obtaining point-wise 
estimates. For further details on (continuous) maximum principles refer to [40 1 115 1 ITH ] 117 1 I^ H I20j . 

A direct consequence of maximum principles for diffusion-type is the non-negativity of the so- 
lution (under appropriate conditions on the source and boundary conditions). Physical quantities 
like concentration of a diffusant should be non-negative by their nature and their approximations 
should also be non-negative as well. The question to ask is whether a chosen numerical formulation 
satisfies these maximum principles and meet the non- negative constraint. The discrete version of 
maximum principles is commonly referred as discrete maximum principles. 

Many existing numerical formulations and packages do not satisfy the maximum principles. They 
may produce negative values for the primary variables in diffusion-type equations (that is, negative 
values for the concentration and temperature). It should be emphasized that the violation is not 
mere numerical noise, that is, the violations will be much larger than the machine precision and 
cannot be neglected. 

For example, in Figure [TJ we have shown that the contours of concentration obtained using 
Abaqus [T] for pure diffusion. The uncoupled problem is similar to one considered in Section 
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HI As one can see from the figure, significant part of the domain has negative solution. The 
minimum value of the concentration is approximately —0.0832, which is 4.16% off the range of 
possible values (which is between and 2). In a subsequent section we show that the classical 
single-field Galerkin formulation produces negative values for the concentration for the same test 
problem (see Figure [15]) . and produces (qualitatively and quantitatively) wrong results for a cou- 
pled deformation-diffusion problem. Furthermore, Nakshatrala and Valocchi [35] have shown that 
the lowest order Raviart-Thomas [48] and variational multiscale [36] H3] mixed formulations violate 
discrete maximum principle and produce negative solutions for pure diffusion equation. Nagara- 
jan and Nakshatrala [52] have shown that the single-field formulation violates discrete maximum 
principles for diffusion with decay. 

1.2. Our approach. Herein, we will consider realistic problems, allow the medium to be inhomo- 
geneous, consider anisotropic diffusion, and develop a robust computational framework that will 
always produce physically meaningful non-negative values for the concentration of the diffusant on 
general computational grids. The computational framework will consist of a non- negative formula- 
tion for diffusion equation, a single-field formulation for the deformation problem, and a staggered 
coupling algorithm. 

We employ a staggered coupling technique (also known as partitioned solution approach) to cou- 
ple individual analyses to obtain the coupled response. The non-negative formulation for diffusion 
is developed by extending the Galerkin formulation using convex programming, which produces 
(physical) non-negative solutions for the concentration even on general computational grids with 
low-order finite elements (linear three-node triangular, bilinear four-node quadrilateral, linear four- 
node tetrahedron, and tri- linear eight-node brick elements). The proposed non- negative formulation 
being applicable only to low-order finite elements is not a limitation as low-order finite elements 
remain quite popular in the solution of practical problems. (This is particularly true for large-scale 
simulations with complex geometries because of the inherent simplicity of low-order elements. An- 
other reason is that adaptive mesh-generation techniques are simpler and tend to perform better 
with low-order finite elements.) 

1.3. Main contributions of this paper. Some of the main contributions of this paper are as 
follows: 

(1) We presented a mathematical model for diffusion of an inert chemical species in an elastic 
deformable solid that takes into account the effect of diffusant on the deformation. The model 
is truly coupled in the sense that the deformation of the solid will be affected by the diffusion 
process, and the diffusion process is in turn affected by the deformation of the solid. Such a 

3 



coupled deformation-diffusion model is suitable to study degradation/healing in elastic solids. 
We restricted our model to steady-state. 

(2) We presented a robust computational framework for performing deformation-diffusion analysis. 
The framework includes a solver for deformation, a non-negative solver for tensorial-diffusion, 
and a coupling algorithm to couple the individual deformation and diffusion analyses. The 
coupling algorithm is a staggered algorithm in which deformation and diffusion sub-problems are 
solved in an iterative fashion until convergence. An important aspect of the proposed framework 
is that it employs a novel numerical formulation that ensures non-negative solutions for the 
concentration of the diffusant. 

(3) Using the proposed computational framework, we solved some realistic finite domain problems 
on general computational meshes. Also, we systematically studied the effect of the concentra- 
tion of the diffusant on the deformation of the solid and vice- versa, and their implications on 
degradation/healing of materials and structures. 



1.4. An outline of this paper and symbolic notation. The remainder of this paper is orga- 
nized as follows. Section [2] presents a mathematical model for degradation/healing of a deformable 
elastic solid. The mathematical model will require performing coupled deformation-diffusion analy- 
sis. In Section[3]we present a fully coupled computational framework for deformation-diffusion anal- 
ysis. The proposed computational framework will contain a non-negative formulation for tensorial- 
diffusion equation on general computational grids (which will always produce physically meaningful 
non-negative values for the concentration), a numerical solver for deformation, and a staggered cou- 
pling algorithm for coupling individual diffusion and deformation numerical solvers. In Section [H 
representative numerical examples will be presented to illustrate the performance of the proposed 
coupled deformation-diffusion computational framework. Conclusions are drawn in Section [SJ 

The symbolic notation adopted in this paper is as follows. We shall make a distinction between 
vectors in the continuum and finite element settings. Similarly, we make a distinction between 
second-order tensors in the continuum setting versus matrices in the context of the finite element 
method. The continuum vectors are denoted by lower case boldface normal letters, and the second- 
order tensors will be denoted using upper case boldface normal letters (for example, vector u and 
second-order tensor T). In the finite element context, we shall denote the vectors using lower 
case boldface italic letters, and the matrices are denoted using upper case boldface italic letters 
(for example, vector v and matrix K). Other notational conventions adopted in this paper are 
introduced as needed. 
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2. A MATHEMATICAL MODEL FOR COUPLED DEFORMATION-DIFFUSION 



Consider an inert (chemical or biological) species being diffused through a deformable elastic 
solid. We now present a simple mathematical model for such a process. Let SI C R nd be a bounded 
open domain, where "racf ' denotes the number of spatial dimensions. The boundary d£l is assumed 
to be piecewise smooth. Mathematically, d£l = — 0, where is the set closure of A spatial 
point in is denoted by x. The gradient and divergence operators with respect to x are denoted by 
grad[-] and div[-], respectively. The unit outward normal to the boundary is denoted by n(x). We 
shall denote the displacement of the solid by u(x), and the concentration of the diffusant by c(x). 
It is important to note that the concentration is a non-negative quantity, and a robust numerical 
solver should meet the non-negative constraint (which is not the case with many popular numerical 
schemes). 

For the deformation problem, the boundary is divided into two complementary parts: on 
which the displacement vector is prescribed, and on which the traction vector is prescribed. For 
the diffusion problem, the boundary is divided into rj? on which concentration is prescribed, and 
on which the flux is prescribed. For well-posedness, we require that rj^nr^ = 0, rj^ur^ = dQ, 
fl rf = 0, and rj? U = <90. In addition, we assume that meas(r^) > and meas(r^) > 
for uniqueness of the solution. 

2.1. A model for diffusion-dependent deformation. The solid is modeled using linearized 
elasticity but the material parameters are allowed to depend on the concentration. The linearized 
strain is defined through 

E; := - (grad[u] + grad[u] T ) (1) 

For a given concentration c(x), the stress-strain relationship will be modeled as follows: 

T c (u, x) = A(x, c)tr[E,]I + 2/x(x, c)E, (2) 

where T c is the Cauchy stress, and A and /x are the Lame parameters but now can depend both on 
the concentration and position. A simple model for the Lame parameters to account for degrada- 
tion/healing of the material due to the presence of a diffusant can be taken as follows: 



A(x,c) = A (x) + A!(x)^ (3a) 

Cref 




Cref 



where Cr e f is the reference concentration (which depends on the problem); Ao and no are the 
Lame parameters for the virgin material (that is, in the absence of the diffusant); and Ai(x) 
and /^i(x) are the weights to account for the effect of concentration on the Lame parameters. 
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The material parameters Ai and \i\ can be individually positive (which means that the material 
is healing), negative (which means that the material is degrading), or zero (which means that 
the material is unaffected by the presence of the diffusant). However, it is assumed that the 
parameters (Ao(x), Ai(x), /^o(x) and //i(x)) and concentration are such that we have bulk modulus 
A(x, c) + |/u(x, c) > and shear modulus /i(x, c) > 0. 

The above model given by equation ([3]) is similar to the concept of macroscopic damage variables, 
which has been introduced by Kachanov [21] to model material damage. The basic idea behind 
the concept of macroscopic damage variables is to quantify damage using internal variable(s). This 
concept has now been widely employed in numerous other works on damage (for example, see the 
review articles by Chaboche [3 [9]). In the above model, the concentration of the diffusant can 
be thought as a macroscopic damage/healing variable. The most common criticism about using 
internal variables is that these variables (in many cases) cannot be measured using physical experi- 
ments. However, in the above model, the dependence of the Lame parameters on the concentration 
can be measured indirectly by non-destructive testing methods. 

It should be noted that Weitsman |56j . Kringos et al. |29} 128] . and Muliana et al. [41] have 
considered material degradation as a function of concentration of the diffusant. However, the 
models developed in these works are not fully coupled. That is, the material properties of the solid 
are dependent on the concentration but the diffusivity is unaffected by the deformation of the solid. 
In Reference [25], a fully coupled model is proposed but a specific boundary value problem (torsion 
of a cylindrical annulus undergoing degradation) is solved using the semi-inverse method. They did 
not present a computational framework, and also their mathematical model is different from the 
one present in this paper (see Remark 12. ip . 

2.2. A model for deformation-dependent diffusion. We define the first and second invariants 
of the tensor E; as follows: 

I El := tr[E,] (4a) 



II El := V2dev[E z ].dev[E z ] = J- (3tr[Ef] - (tr[E,]) 2 ) (4b) 



where dev[E;] := E; — |tr[E/]I is the deviatoric part of Ej. (Note that the second invariant is not 
the principal invariant, and the reason for such a choice will be discussed later.) For a given strain 
field (that is, for a given deformation field), we model the effect of deformation on the diffusivity 
as follows: 

D Ei (x) = D (x) + (D T (x) - D (x)) ( ^^4) + (D s (x) - D (x)) ( ^MlMzl) 

V ex P far-Eref ] - 1 / V ex P [r]sE re{ ] - 1 / 

(5) 



where t]t and rjs are non-negative parameters; Do(x), Dt(x) and Dg(x) are (respectively) the 
reference diffusivity tensors under no, tensile, and shear strains; and E rc f is a reference measure 
of the strain. The above model is partly motivated by the stress-induced diffusion experiments on 
glass done by McAfee |38l 139]. These experiments have clearly shown the following aspects, which 
have been qualitatively incorporated in the above model. 

(a) The relative diffusion rate under tension is nearly five times more than that of the relative 
diffusion rates under compression and shear. 

(b) The relative diffusion rate varies exponentially with respect to the (circumferential) strain for 
Pyrex glass (see [38, Equation 15 and Figure 3]). 

(c) The relative diffusion rate under compression is significantly different from that of shear (see 
[391 Figure 4]). 

Remark 2.1. In Reference [25], the effect of deformation on the diffusivity tensor is modeled using 
the Frobenius norm of the Almansi-Hamel strain. (In linearized elasticity, the Almansi-Hamel 
strain is approximately equal to linearized strain.) Although the Frobenius norm of the strain is an 
invariant, the model cannot capture the difference in diffusivity under tension, compression, and 
shear, which has been observed in many materials. On the other hand, our model given by equation 
(0) can capture such departures between the diffusivities under tension, compression and shear. 

Remark 2.2. A remark on the choice of invariants in the model (given by equation §5$)) is in 
order. Note that the second principal invariant of a tensor A is defined as 



which is different from the one used in equation @ . It has been discussed in the literature that 
the principal invariants are not suitable to fit experimental data (e.g., Lurie [33], Anand [21 [3J, 
Criscione et al. [13], Plesek and Kruisovd j49j ). These works employed invariants of Eulerian 
Hencky strain. Let the Eulerian Hencky strain be denoted by E# = ln[V], with V being the left 
stretch tensor in the polar decomposition of the deformation gradient F. The first three invariants 
based on the Eulerian Hencky strain which represent dilation (k\), magnitude of distortion fa) and 
mode of distortion fa ) are given by 




(6) 



fci 



tr[E#] = ln[J] 





(7b) 



3\/6det — dev[En] 




'2 
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where dev[E#] := E# — ^tr[E#]I and J = det[F] > 0. For small gradients of the displacement, 
we have 

Eh ~ E; (8) 

In linearized elasticity, the constitutive equation depends only on k\ and &2 as any dependence 
on &3 makes the model inherently non-linear (see Criscione et al. [13] )■ Herein, for modeling 
deformation- dependent diffusivity tensor we did not use k^ just to be consistent with the theory of 
linearized elasticity. However, in the case of finite elasticity, one should also use the third invariant 
in modeling the deformation- dependent diffusivity. 

One of the attractive features of the above model is that the parameters in the model can be 
established using standard experiments, which we shall describe below. 

(a) It is easy to check that if there is no strain (that is, E; = 0), then De, = Do- Hence, one 
can find the reference diffusivity tensor Do by doing a diffusion experiment on an unstrained 
specimen. 

(b) Under simple shear, we have Je; = 0, and IIe 1 = as, where as is the angle of shear. As 
as — >■ E rc f, the diffusivity tensor D — > D5. Hence, one can find Dg and r/s by doing a diffusion 
experiment on a specimen undergoing a simple shear. 

(c) Under a uniform tri-axial tension test, the linearized strain can be written as E^ = a^I with 
ar > 0. In this case, 11^ = 0, Ie ; = 3ot, and as aj — > ^E Te { we have D — > T>t- Hence, one 
can find Dy and r/T in the model given by equation Q by doing a diffusion experiment on a 
specimen under uniform tri-axial strain. 

There are experimental techniques discussed in the literature for maintaining a specimen 
under uniform tri-axial strain. To name a few, tri-axial tensile testing of brittle materials 
such as calestone (a dental plaster), copper and aluminum alloys, austenitic stainless steel were 
described by Cridland and Wood [12], Hayhurst and Felce [21] and Calloch and Marquis [7]. 
Advanced tri-axial testing of geomaterials such as rock and soil were carried out by Donaghe 
et al. |14j . Hunsche [23] and Wawersik [S3]. Despite these experimental techniques, it can well 
be argued that maintaining a specimen under uniform tri-axial strain can be a difficult and 
expensive. In that case, after determining Do, D5 and rjs one can evaluate T>t and rjr using 
either uni-axial or bi-axial tension tests. However, in uni-axial and bi-axial tests, it should be 
noted that J/e ; will not be equal to zero. 

Remark 2.3. For brittle materials such as concrete, ceramics, metallic alloys and geomaterials 
such as soil and rock it is relatively easier to perform compression tests than tension tests. In those 
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cases, the expression for diffusivity tensor D E; (x) can be modeled as follows: 
D El (x) = D„(x) + (D (x) - D C (X)) ( ^E,]-l \ + _ / expfejft ] - A 

(9) 

where rjc an d are non-negative parameters; Do(x), Dc(x), and Dj(x) are (respectively) ref- 
erence diffusivity tensors under no, compressive, and shear strains. Also, it should be noted that 
modeling diffusivity as an exponential function of stress/strain is quite popular in literature (e.g., 
see McAfee [38], Fahmy and Hurt |16| ). 

The diffusivity tensor is assumed to be symmetric, bounded, and uniformly elliptic. That is, 

D Ei (x) =D| ( (x) VxgO (10) 
and there exists two constants < £i < £2 < +00 such that 

£iy T y < y T D E; (x)y < £ 2 y T y VxeH and V y e R nd (11) 

2.3. Governing field equations. The governing equations for the deformation of the solid can 
be written as follows: 

div[T c ] + p(x)b(x) =0 in O (12a) 
u(x) = u p (x) onT° (12b) 
T c n(x)=tP(x) ouT* (12c) 

where p(x) is the density, b(x) is the specific body force, u p (x) is the prescribed displacement, 
t p (x) is the prescribed traction, and recall that n(x) is the unit outward normal to the boundary. 
In the absence of internal couples, the balance of angular momentum reads 

T c = T T C (13) 

which the Cauchy stress given in equation (J2J clearly satisfies. The governing equations for the 
steady-state (deformation-dependent) diffusion process can be written as follows: 

-div [D E; (x) grad[c]] = /(x) in (14a) 
c(x) = c p (x) on T° (14b) 
n(x) • D Ei (x) grad[c] = /i p (x) on if (14c) 

where c p (x) is the prescribed concentration, /i p (x) is the prescribed concentration flux on the 
boundary, and /(x) is the volumetric source. 

It is easy to see that the governing equations for the deformation (|12p and the governing equations 
for the diffusion (fT4"|) are coupled through equations and To predict the degradation/healing 



of the solid due to the diffusant, one needs to solve this system of coupled nonlinear partial dif- 
ferential equations given by equations (0), (fT2|) and {HD- It is noteworthy that, except for 
simple problems, it is not possible to find analytical solutions to this system of equations, and one 
needs to resort to numerical solutions for solving realistic and practical problems. 

However, to obtain reliable and predictive numerical solutions, one has to overcome many nu- 
merical challenges. In particular, one has to make sure that the chosen numerical scheme gives 
non-negative values for the concentration as a negative value for the concentration is unphysical. 
We will show in a subsequent section that the classical single-field formulation produces negative 
concentrations. This is particularly true if the medium is anisotropic, and one may even get nega- 
tive values for the concentration even when the medium is isotropic if the mesh is not chosen with 
care. For example, one need to choose a mesh with square elements or a well-centered triangular 
mesh [45| even for isotropic diffusion. (In two-dimensions a well-centered triangular mesh means 
that all the angles of any triangle are acute. Similarly, one can define a well-centered mesh in higher 
dimensions [53].) 

We now present a numerical framework to solve the coupled deformation- diffusion equations (fTJ) 
(|14p in a systematic manner. 



The computational framework will be based on the Finite Element Method (FEM), convex 
quadratic programming, and staggered coupling techniques. To this end, let the domain £1 be 
decomposed into "iVeZe" non-overlapping open sub-domains (which in the finite element context 
will be elements). That is, 



where a superposed bar denotes the set closure. The boundary of Q e is denoted as d£l e := £l e — £l e . 
For convenience and to avoid errors due to projection operators, we shall employ the same com- 
putational mesh for both deformation and diffusion analyses. (Note that one needs to employ 
projection operators if different computational meshes are employed for multi-field problems like 
coupled deformation-diffusion.) We now present individual solvers for deformation and diffusion, 
and a coupling algorithm to couple individual solvers to obtain the coupled response. The solver 
for the diffusion problem will always give physically meaningful non-negative values for the con- 
centration. 



3. A COUPLED COMPUTATIONAL FRAMEWORK 



Nele 




(15) 



e=l 
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3.1. A numerical solver for deformation analysis. For solving the pure deformation problem 
(that is, for a given concentration field) the computational framework utilizes the standard single- 
field (pure displacement) formulation. However, it should be noted that one can employ any other 
formulation (e.g., B-bar method [52], stabilized mixed formulation [33]; mixed assumed strain 
formulation [50], mixed enhanced formulation [52] ) to solve the (pure) deformation problem. For 
completeness, we shall now briefly outline the single-field formulation. To this end, we shall define 
the following function spaces: 

U := ju(x) G (H 1 (Q)) nd | u(x) = u p (x) on T°} (16a) 

W := {w(x) G (H 1 (n)) nd | w(x) = on T°} (16b) 

where H 1 ($7) is a standard Sobolev space on O 0, and recall that "nd" is the number of spatial 
dimensions. The standard single-field formulation for the pure deformation problem (|12|) reads: 
Find u(x) G U such that we have 

£ u (w;u) = L u (w) Vw(x)GW (17) 

where the bilinear form and linear functional are, respectively, defined as follows: 

B u (w;u):= / gradfw] • T c (u,x) dtt (18a) 
Jn 

L u (w) := [ w(x) • p(x)b(x) dft + f w(x) • t p (x) dr (18b) 
Jn Jv* 

For a non-negative integer m, let P m (f2 e ) denotes the linear vector space spanned by polynomials 

up-to mth order defined on the sub-domain J7 e . We shall define the following finite dimensional 

subsets of U and W: 

U h := |u h (x) G U | u fe (x) G (C°(ti)) nd , u fc (x)| ne G (V 1 ^ 6 ))™*, e = 1, • • • ,iVe/e| (19a) 

W h := |w fe (x) G W | w^x) G (C°(n)) nd , w /l (x)| ne G (ff> fe (ft e )) ^ , e = 1, • • • ,iVe/e| (19b) 

where k is a non-negative integer. A corresponding finite element formulation for the deformation 
analysis can be written as: Find u' l (x) G U h such that we have 

B u {w h ; u h ) = L u {w h ) V w /l (x) G W l (20) 

After the finite element discretization, the deformation analysis will involve solving the following 
discrete equations: 

K u {c)u = f u (21) 

where u denotes the nodal displacements, c denotes the nodal concentrations, and the stiffness 

matrix for the deformation analysis K u depends on the concentration of the diffusant. 

li 



3.2. A numerical solver for diffusion analysis. Before we provide a numerical solver for dif- 
fusion, we shall provide a mathematical argument to show that c(x) > in Cl even for the coupled 
problem. We shall assume that a solution exists for the coupled problem. (Proving existence of a 
solution for the coupled system is a research topic by itself, and is beyond the scope of this paper. 
However, in this paper we do find a numerical solution for various coupled deformation-diffusion 
problems.) That is, there exists a pair, u(x) and c(x), such that they satisfy the coupled system 
of equations. For the given displacement field, u(x), (and hence for a given strain field E/(x)) we 
define 



From the theory of partial differential equations we have the following maximum principle [15j : 

Theorem 3.1. Let c p (x) > on dCl and D(x) be continuously differentiable. If c(x) G C 2 (Cl) n 
C°(Cl) satisfies the differential inequality — div[D(x)grad[c]] = /(x) > in Cl, then we have the 
following non-negative property: 



Remark 3.2. The above maximum principle theorem is due to Hopf, and a proof can be found in 
any standard textbook on partial differential equations (e.g., References |4C H 115 1 I^ H H8| . I20j ). One 
can find in the literature maximum principles for diffusion-type equations under weaker regularity 
than C 2 (Cl) n C° '(CI) (for example, see References [371 H]J- But such a thorough treatment is beyond 
the scope of this paper, and is not crucial to our presentation. 

As discussed earlier, many existing numerical formulations (including the single-field formulation, 
which is based on the Galerkin principle) for diffusion-type equation do not meet the non-negative 
constraint. For example, the widely used single-field formulation (which is based on the Galerkin 
principle) does not produce physically meaningful non-negative solutions. We now present a novel 
methodology of enforcing the non- negative constraint on the concentration of the diffusant. The 
methodology works well for general computational grids and for low-order finite elements. 

We start with the single-field formulation, and then modify the underlying (discrete) variational 
statement to meet the non-negative constraint. We shall define the following function spaces 




(22) 




(23) 




(24a) 



Q:={q(x)eH 1 (Cl)\q(x)=OonT^} 



(24b) 
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where H l (tt) is a standard Sobolev space [6]. We also relax the regularity of the diffusivity tensor 
for weak solutions, and assume that 

/ tr[D Ej (x) r D Ei (x)] dft < +00 (25) 
Jn 

where tr[-] is the standard trace operator used in tensor algebra and continuum mechanics |10j . 
The standard single field formulation for tensorial diffusion equation (|14p reads: Find c(x) G V 
such that we have 

B c (q;c)=L c (q) V g(x) G Q (26) 

where the bilinear form and linear functional are, respectively, defined as 

B c (q;c):= [ gmd[q] ■ D E; (x)grad[c] dft (27a) 
Jn 

L c (q) := [ q(x) /(x) dft + f ? ( x ) fcP( x ) dr (27b) 

It is well-known that the above weak form (|26p has a corresponding variational statement, which 
can be written as follows: 

minimize -B c (c;c) — LJc) (28) 
c(x)eP 2 

We shall define the following finite dimensional vector spaces of V and Q: 

V h := |>(x) G P I c h (x) G C°(fi), c ft (x)| ne G P fe (ft e ), e = 1, • • • , iVe/e} (29a) 
Q h := |V(x) G G I ^(x) G C°(fi), g h (x)| ne G F k (Q e ), e = 1, ■ ■ ■ ,Nele} (29b) 

where A; is a non- negative integer. (Recall that, for a non-negative integer m, P m (O e ) denotes the 
linear vector space spanned by polynomials up-to mth order defined on the subdomain il e .) A 
corresponding finite element formulation can be written as: Find c ft (x) G V h such that we have 

B c (q h -c h ) = L c (q h ) V/(x)eQ k (30) 

3.2.1. A non-negative solver for tensorial diffusion equation. We shall use the symbols ^ and >z. 
to denote component-wise inequalities for vectors. That is, for given any two (finite dimensional) 
vectors a and b 

a <b means that dj < 6j V i (31) 

Similarly, one can define the symbol y. We shall denote the standard inner-product on Euclidean 
spaces by (•; •). 
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After finite element discretization, for given nodal displacement vector u, the discrete equation 
for the diffusion analysis takes the following form: 



K c {u)c = f c (32) 

where K c is a symmetric positive definite matrix, c is the vector containing nodal concentrations, 
and f c is the source vector. Let "ndofs" denote the number of degrees-of-freedom for the concen- 
tration. The matrix K c is of size ndofs x ndofs, and the vectors c and f c are of size ndofs x 1. 
The finite element discretized equation (|32p is equivalent to the following minimization problem: 

minimize -(c;K c (u)c) - (c; f c ) (33) 



As shown in Figure 15(a) , the numerical formulation based on equations f)32 [) and (|33p produces un- 
physical negative concentrations for many practically important problems. (More examples showing 
Galerkin formulation producing negative solutions can be found in Reference |42j.) Following the 
ideas outlined in References |45} l32j 02] a non-negative formulation corresponding to (J33]) can be 
written as follows: 

minimize — (c; K c (u)c) — (c; f c ) (34a) 



subject to c h (34b) 

Since for a given nodal displacement vector, the matrix K c (u) is positive definite, the above problem 
(|34j) belongs to convex quadratic programming. From optimization theory it can be shown that 
the problem (|34p has a unique global minimizer. There are many robust numerical algorithms 
available in the literature to solve the aforementioned constrained minimization problem (e.g., 
active set strategy, interior point methods, barrier methods). In all our numerical simulations, we 
have employed the active set strategy. A detailed discussion on numerical optimization can be 
found in references 061 [M] . 



3.3. A staggered coupling algorithm. Current coupling algorithms are broadly classified into 
main classes: monolithic and staggered schemes. In monolithic schemes discretization scheme is 
applied to the full problem. On the other hand, in staggered schemes (which are based on operator- 
split techniques) the coupled system is partitioned, often according to the different coupled fields 
(in our case, concentration and displacement in deformation-diffusion analysis), and each partition 
is treated by different and tailored numerical schemes. There is no easy way to incorporate our non- 
negative formulation within a monolithic scheme. Therefore, we shall employ a staggered coupling 
approach. The various steps of the proposed coupling algorithm is presented in Algorithm [TJ 
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Algorithm 1 Staggered coupling algorithm for deformation-diffusion analysis 



1: Input: tolerance (&j?ol)' max i mum number of iteration (MAXITERS) 

2: Guess c(°) t 

3: for i = 1, 2, ■ ■ ■ do 

4: if t > MAXITERS then 

5: Solution did not converge in specified maximum number of iterations. EXIT 

6: end if 

7: Call deformation solver: Obtain «W by solving 

8: Call non-negative diffusion solver: Obtain C" by solving the following minimization problem 

minimize -(c®; Jjf c («W)cW> - (c^;/ c ) 

c (i)gK™ d °/ s 2 

subject to c w ^ 



9: if ||cW - c^" 1 ))! < e^' QL then 
10: Staggered scheme converged. EXIT 
11: end if 
12: end for 



Remark 3.3. In this paper, we shall use the 2-norm in the stopping criterion ||c® — c^ _1 ^|| < 
e TOL ^ n staggered coupling algorithm. However, one could use any other norm as in (finite 
dimensional) Euclidean spaces all norms are equivalent |19j . 

Remark 3.4. It should be noted that even though the solution procedure is a staggered scheme, the 
problem is coupled, and the converged numerical solution will be the coupled response. 

4. PERFORMANCE OF THE PROPOSED FRAMEWORK 

In this section, we present numerical solutions of several realistic problems using the proposed 
mathematical model and computational framework. The first set of problems involve degrada- 
tion/healing of beams, and the second set involves degradation/healing of rectangular domains 
with holes. All these problems naturally arise in many important engineering applications. For 
example, study of degradation of beams is important to assess the performance of structural com- 
ponents in bridges, towers, buildings and aircraft wings. The second set of problems is to study 

the degradation of ducts carrying chemicals (e.g., water, carbon-dioxide, coolant) cutting through 
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a slab or a wall, which have many applications in civil and nuclear structures. Using these two 
sets of representative problems we will illustrate the performance of the proposed computational 
framework for coupled deformation-diffusion analyses. In particular, we will show that the proposed 
computational framework produces physical and reliable solutions. We also systematically study 
the effect of the concentration of the diffusant on the deformation of the solid and vice-versa. 
Here, in these representative numerical examples we have chosen Do, Dy and D5 as follows: 

Dq= / cos(fl) -sin(0) W cos(0) sin(0) \ 



sin(0) cos(0) / V d 2 J \ -sin(0) cos(6») 
D T = <£ T D (35b) 
D 5 = $ 5 D (35c) 

where $t and &s are some positive real numbers. The reason for choosing such a form is that 
a change in each of the Dy and D5 effect the concentration in a considerable manner in various 
important realistic problems. But in general each component of Dj-, D5 may be different from that 
of Do, even then our computational framework is still valid and works as shown in the representative 
problems outline in subsections plate with a hole 14.5.11 and beam with three holes I4.5.2| where in 
D7 1 and Dg are chosen as follows: 

/ cos(0) -sin(0) W df W cos(0) sin(0) \ 
I sin(0) cos(#) J \ dj / V - sin(0) cos(#) / 




/ cos{9) _ sin(e) \ / d s o W cos{6) s[n{9) \ 
S \ sm(6) cos(#) J \ df / \ -sin(0) cos(6) J 

4.1. Numerical /i-convergence study. In this subsection, we perform a numerical /i-convergence 
study by employing the proposed computational framework on a coupled deformation-diffusion 
problem. Herein, by /i-convergence we mean the overall convergence of the proposed framework 
with respect to the refinement of the computational mesh (but still maintaining the same order of 
interpolation within each finite element). Since the non- negative solver for diffusion works only for 
low-order finite element, we employ low-order finite elements in all our numerical simulations. It 
should be noted that for each successful coupled analysis using the proposed computational frame- 
work, the staggered coupling algorithm should converge. For each iteration in the staggered coupling 
algorithm, the active-set strategy in the non-negative solver for the diffusion problem should con- 
verge. 

We employ the method of manufactured solutions |27] in this convergence study. The computa- 
tional domain is taken as a bi-unit square (0 < x < 1 and < y < 1). The displacement vector is 
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given by 



The concentration is given by 



u{x,y) = -sm(— )sm( — ) 

7T 2 Z 

1 7tx ny 
v(x,y) = -cos(— )cos(— ) 

7T 2 I 



c{x, y) = 1 + - sm( — ) sin(-) 

7T 2 2 



(37a) 
(37b) 

(38) 



Note that the concentration given in the above equation is non-negative in the whole computational 
domain. The following parameters are assumed in the convergence study: 



Ho = vr + 2, ^ = -7r, A = 2, Ai = 

Cref = 1, f]T = 1, VS = 1, Pt = 2, l3 S = 2, 7s 

D = 2 I, D T = /3 T D , D 5 = /3 5 D 



1, P = 1, ^rcf = 0.0001, 

Ps-l 



exp[rj s E ief ] - 1' 



(39) 



The volumetric source for the diffusion problem is given by 



f{x,y) = vrsin(— )sin(- 



(1 - is) +75 exp 



irx . ny i\ 
cos(— )sm(— )jj 



— — sm(7rx) cos(7r?/jexp cos(— )sm(— ) 



(40) 



The specific body force for the deformation problem is given by 
b(x,y) = 



(41) 



(7r sin(^) sin(^) + j cos(7rx)(l — cos(7ry)) \ 
7r cos( 2 y ) cos(^) — j sin(7ra;) sin(7ry) J 

The boundary value problem is pictorially described in Figure [2j The tolerance in the stopping 
criterion for the staggered coupling algorithm is taken as e^OL = 10 -8 - A hierarchy of meshes 
similar to the ones shown in Figure [3] is employed in the numerical convergence study. In Figure 



HI the obtained numerical solution using the mesh shown in Figure 3(a) is compared with the 



analytical solution. The convergence of the proposed coupled framework is illustrated in Figures [5] 
and [6] for three-node triangular and four-node quadrilateral elements with respect to L 2 -norm and 
ff 1 -seminorm, and the proposed computational framework performed well. 

In the subsequent subsections, we employ the proposed computational framework to solve various 
finite domain practical problems. Using these we illustrate how the tension, compression and shear 
in the solid affect the (steady-state) diffusion process, and how the concentration of the diffusant 
affect the deformation of the solid. The values of the parameters that are common to all the test 
problems are presented in Table [TJ We have given the values for the volumetric source and body 
force for all the test problems in Table [2J 
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Table 1. Values of the parameters that are the same in all test problems (except for the 
one in the numerical /i-convergence). 



Parameter 


value 


P 


i 
± 


Cref 


1 


^ref 


0.0001 


Ao 


+10 6 


/'o 


+10 6 


Ai 


-9 x 10 5 


Mi 


-9 x 10 5 



Table 2. Volumetric source /(x) and body force b(x) for various test problems. 



Test problem /( x ) b(x) 



Cantilever beam with edge shear 


10000 







Simply supported beam under self-weight 


1000 


-10 


e» 


Fixed beam under self- weight 


100 


-10 




Plate with a hole under self-weight 





-10 




Beam with three holes under self-weight 





-10 


e y 



4.2. Cantilever beam with edge shear. The purpose of this test problem is to illustrate the 
affect of Ds on the coupled response. We consider a cantilever beam with a uniform edge shear 
of 500. The length of the beam is 1.0, and the height of the beam is 0.1. For the deformation 
problem, top and bottom surfaces are subjected to zero traction. For the diffusion problem, zero 
concentration is prescribed on the top and bottom surfaces, and zero flux is prescribed on the left 
and right surfaces. The boundary value problem is pictorially described in Figure 

We have employed a structured mesh using four-node quadrilateral elements with 21 nodes 
along each side of the domain (see the mesh in Figure [8j). The stopping criterion is again taken 
as e^oL = 10 -8 . We have considered three different values for <£,$: 5, 10, and 20. The other 
parameters used in the coupled analysis are as follows: 

9 = J, dx = 1, d 2 = 1, $ T = 2, VT = 100, r]s = 1 (42) 
o 

The contour profiles for the concentration for these three values are shown in FigureO As expected, 
the diffusant accumulated near the right side of the domain (where the uniform edge shear is 
applied). In Tabled we have presented the maximum concentration and the number of iterations 
taken by the staggered coupling algorithm for the chosen values of <J>s. One can observe that as 
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Table 3. Cantilever beam with edge shear: The table presents the maximum concentration 
and the number of iterations taken by the staggered coupling algorithm for various values 
of $5. (Note that the minimum concentration in all the cases is zero.) 



$5 Max. concentration Iterations 
5 4.257 x KT 1 14 
10 2.187 x 10- 1 9 
20 1.107 x 10" 1 7 



Table 4. Simply supported beam under self- weight: The table presents the maximum 
concentration and the iterations taken by the staggered coupling algorithm for various values 
of rjs- (Note that the minimum concentration in all the cases is zero.) 

rjs Max. concentration Iterations 
1 7.205 x 10" 1 10 
10 3 7.309 x 10- 1 10 
2 x 10 4 9.365 x 10" 1 12 



&s increases, the staggered coupling algorithm takes fewer iterations. This can be explained as $5 
increases, the diffusivity increases and the diffusant is more uniformly distributed (which is evident 
from Figure [8J . 

4.3. Simply supported beam under self-weight. The purpose of this test problem is to study 
the affect of r]s on the coupled response, and to illustrate the competitive effects of shear and 
tension/compression on the diffusion process. We consider a simply supported beam under self- 
weight with traction-free surfaces. The boundary value problem is pictorially described in Figure 

m 

(c) 

The tolerance for the stopping criterion in the staggered coupling algorithm is taken as c^ql = 
10~ 5 . We have employed a structured mesh using four-node quadrilateral element with 21 nodes 
along each side of the domain (see the mesh in Figure [TUJ). We have considered three values for rjs- 
1, 10 3 and 2 x 10 4 . The other parameters used in this test problem are as follows: 

9 = J, di = 1, d 2 = 1, $ T = 10, $ s = 10, T) T = 1 (43) 

From Table U] one can infer that as r]s increases the maximum concentration and iterations taken 

for convergence of staggered coupling algorithm increase. This is because for lower values of rjs 

shear affects the diffusivity tensor more than the tension/compression and for higher values of rjs 

tension/compression dominates. This is also evident in the numerical results presented in Figure 
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Table 5. Fixed beam under self- weight: The table presents the maximum concentration 
and the iterations taken by the staggered coupling algorithm for various values of T]t- (Note 
that the minimum concentration in all the cases is zero.) 



<3?T 


Max. concentration 


Iterations 


1 


1.250 x KT 1 


2 


5 


1.348 x KT 1 


5 


7 


1.575 x KT 1 


8 



[PHI For lower values of r]g, the diffusant accumulated near the supports of the beam (where the 
shear is maximum). For higher values of r/s, the diffusant spreads deep into beam. An important 
feature to be noted is that the concentration profile curves up. This is because the diffusivity is 
higher in tension than in compression. This is physically meaningful, as, in general, the sizes of 
pores in the solid enlarge due to tension, shrink due to compression, and distort due to shear. 
The deformation-dependent diffusivity tensor De, (x) given by equation ([5]) takes these factors into 
account. The next test problem highlights other important features of the proposed mathematical 
model. 

4.4. Fixed beam under self- weight. The purpose of this test problem is to study the affect of 
<3?t on the coupled response and show that our model can capture the qualitative aspects of the 
experiments done by McAfee [38} I39j . A fixed beam of length unity and depth 0.1 is subjected 
to self-weight. The top and bottom surfaces are traction-free. For the diffusion problem, the top 
and bottom surfaces have zero concentration, and the left and right surfaces have zero flux. The 
boundary value problem is pictorially described in Figure [TT1 

(c) 

The tolerance for the stopping criterion in the staggered coupling algorithm is taken as e^OL = 
10~ 7 . We have employed a structured mesh using four-node quadrilateral element with 21 nodes 
along each side of the domain (see the mesh in Figure [T2|) . We have considered three different 
values for 1> 5 and 7. The other parameters used in this test problem are as follows: 

9 = 0, di = 1, d 2 = 1, $s = l,rrr = 100, rj s = 1 (44) 

In Table [5j we have presented the maximum concentration and number of iterations taken by 
the staggered coupling algorithm for chosen values of <3?t- One can observe that as $>t increases the 
maximum concentration and iterations taken increase. This results in the deformation-dependent 
diffusivity De ; being high in tensile and low in compressive and shear regions of the beam. The 
diffusant accumulates in low diffusive areas such as the supports which are under shear and top 
section of the beam which is in compression. From Figure [121 one can observe that the concentration 
profile which is initially around the central line of the beam curves up (the region of interest being 
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around 0.12). This test problem shows the success of our mathematical model in predicting the 
underlying physical phenomena observed from the stress-induced experiments done by McAfee 
[38[ [39] . The invariants which represent dilation and distortion in our model were instrumental 
in capturing the change in diffusivity under tension, compression and shear. Hence concentration 
curves towards the compression zone. 

This type of curving of the concentration profile is not observed when one uses the model in which 
De,(x) depends only on ||E/||. As mentioned in the Remark 12, 1\ a diffusivity model based on the 
frobenius norm of E^ cannot capture the change in diffusivity tensor under tension, compression 
and shear deformations. The diffusivity tensor as per Reference [25] and the parameters assumed 
for the numerical study are given by: 

D Ei (x) = D (x) + (D oc (x) - D (x)) (1 - exp[-A||E,||]) (45a) 
= 0, di = 1, d 2 = 1, A = 100, Doo = 10D (45b) 

The concentration profile as illustrated in Figure [13] is always around the central line and does not 
curve up for any value of and A. 

4.5. Degradation/healing of rectangular domains with holes. In the following subsections, 
we highlight the importance of non-negative solutions and its impact on coupled deformation- 
diffusion analyses of rectangular domains with holes. We compare the non-negative formulation 
with the standard Galerkin formulation using two representative problems. We shall model Dy 
and Dg by the equations (j36[) in which d*- is different from dj where i = T, S and j = 1, 2. In both 
numerical studies here, the tolerance for the stopping criterion in the staggered coupling algorithm 
is taken as e^OL = 5 an d three-node triangular unstructured meshes were used (see Figures [T4l 
and [T8l) . 

4.5.1. Plate with a hole under self-weight. The purpose of this test problem is to study the im- 
portance of non-negative constraint in coupled-deformation-diffusion analyses. Herein, we perform 
coupled analysis for a square plate with a square hole under self- weight. The computational domain 
f2 is the region in-between a bi-unit square plate and a square hole of length g. The boundary con- 
ditions for displacements at the hole and traction for the plate are equal to zero. The concentration 
at the hole is maintained at 1 while that at boundary of the plate is equal to 0. The boundary 
value problem is pictorially described in Figure O The other parameters assumed in the analysis 
of the coupled problem el l C cl S follows: 

Q = 1 rjr = 1, ris = 1, di = 10000, d\ = 11000, df = 11000, d 2 = 1, d\ = 10, df = 5 (46) 
3 

From Figure[T5lit is evident that the proposed non-negative formulation satisfies the above condition 
and produces physically meaningful concentration while the standard Galerkin formulation predicts 
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Table 6. Plate with a hole under self- weight: The table presents the minimum concen- 
tration, degradation index and the iterations taken by staggered coupling algorithm using 
the standard single- field formulation. Analysis is done for various values of (di, dj , df ) by 
keeping other parameters in equations 05] fixed. 



(di,cM) 


Min. concentration 


Degradation Index 


Iterations 


Iterations 






(% of nodes violated) 


(Galerkin) 


(non-negative) 


(1,10,10) 








9 


9 


(10,30,20) 








8 


8 


(100,120,110) 


-3.301 x 10~ 3 


25.15 


6 


21 


(1000,1200,1100) 


-3.586 x 10~ 2 


32.76 


5 


22 


(10000,11000,11000) 


-4.398 x 10" 2 


34.07 


4 


8 



otherwise. The comparison of the concentration and degradation profiles for the standard Galerkin 
vs. the non-negative formulation given in Figures [TBI and IT61 shows that the Galerkin formulation 
violates the non- negative constraint. The obtained negative values for the concentration are not 
just numerical noise (see Table [6]). 

4.5.2. Beam with three holes under self-weight. The purpose of this test problem is to show that 
as number of holes increase, the spatial extent and magnitude of violation of the non-negative con- 
straint by the Galerkin formulation increase dramatically. Herein, we perform coupled deformation- 
diffusion analysis for a beam with three square holes under self-weight. The computational domain 
Q is the region in-between a beam of length 10.0 and height 1.0 and three square holes each of 
length 0.4. The boundary conditions for displacements at the holes and traction for the beam are 
equal to zero. The concentration at the holes is maintained at 1 while that at boundary of the 
beam is equal to 0. The boundary value problem is pictorially described in Figure [TTJ The other 
parameters assumed in the analysis of the coupled problem ell C cXS follows: 

e = j t rfr = l,r]s = 1, di = 10000, d\ = 20000, df = 15000, d 2 = 1, d\ = 5, df = 2 (47) 

From Figures [18] and [19] it is evident that as the number of holes increase, the regions of negative 
concentration obtained from the standard Galerkin formulation also increases. Also from Tables [6] 
and [7] it is evident that the negative concentration obtained from numerical study on beam with 
three holes is 1.35 times that of the plate with a hole. 

From Tables [6] and one can see that as the directional diffusivities d\ , df and df increase the 
minimum concentration (which is initially zero) becomes negative under the Galerkin formulation. 
This negative value of concentration and the degradation index (which represents the % of nodes at 
which concentration is negative) increases as these diffusivities increase. The proposed non-negative 
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Table 7. Beam with three holes under self- weight: The table presents the minimum con- 
centration, degradation index and the iterations taken by staggered coupling algorithm using 
the standard single- field formulation. Analysis is done for various values of (di, dj , df ) by 
keeping other parameters in equations 27J fixed. 



(di,cM) 


Min. concentration 


Degradation Index 


Iterations 


Iterations 






(% of nodes violated) 


(Galerkin) 


(non-negative) 


(1,10,10) 








8 


8 


(10,50,30) 








8 


8 


(100,200,150) 


-3.150 x 10~ 2 


30.69 


6 


8 


(1000, 2000, 1500) 


-5.665 x 10~ 2 


31.14 


6 


6 


(10000,20000,15000) 


-5.948 x 10~ 2 


31.14 


6 
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formulation does not give negative values for the concentration for any values for these how these 
diffusivities. This is quite important because the material properties of the solid are dependent on 
the concentration obtained. From Figures [TBI and \19\ one can observe that the standard Galerkin 
method predicts that some regions are healing due to the negative values of the concentration which 
is physically unrealistic. The non-negative solver employed always gives non-negative concentration 
and shows that material is degrading everywhere. 

4.6. Performance of the staggered scheme and the active-set strategy. The convergence 
histories of the staggered coupling algorithm for the aforementioned test problems are shown in 
Figure [20l Note that the convergence of the staggered coupling algorithm is monotonic for the 
test problems on degradation/healing of beams. But the convergence is not monotonic for the test 
problem on degradation/healing of rectangular domain with holes, which is illustrated in Figure 
[20(b) [ ). Figure [2T1 shows the number of active-strategy iterations required for solving the constrained 
optimization problem for each iteration in the staggered coupling algorithm. 

5. CONCLUSIONS 

Many technologically and biologically important processes are coupled deformation-diffusion 
problems. Lately, there is a surge in research activity in studying coupled deformation-diffusion 
problems. However, in all these research efforts it has not been recognized that many popular 
numerical formulations and existing computational packages predict unphysical negative solutions 
for the concentration of the diffusant, and this is more prominent in the case of a medium which 
has high directional diffusivities. Concentration of a chemical is a non- negative quantity, and this 
property has to be preserved to obtain physically meaningful numerical solution for a coupled 
deformation-diffusion problem. 
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We proposed a mathematical model for coupled deformation-diffusion analysis, which is suitable 
to model degradation/healing of materials and structures. The model is fully coupled in the sense 
that the deformation process is affected by the diffusion process, and the diffusion process is in turn 
affected by the deformation of solid. One of the main contributions is that we have presented a 
robust computational framework for solving coupled deformation-diffusion problems. The computa- 
tional framework includes a non-negative formulation for (tensorial) diffusion equation, a numerical 
solver for the deformation of the solid, and a staggered coupling algorithm. An important aspect of 
the computational framework is that it always produces physically meaningful non-negative values 
for the concentration of the diffusant on any computational grid (with low-order finite elements) 
even for a medium which has high directional diffusivities. We have illustrated the robustness of the 
computational framework on representative numerical examples. We have studied systematically 
the affect of the deformation on the diffusion process on various practically important problems. 
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Figure 1. ABAQUS simulation for plate with a square hole: The left figure shows three- 
node triangular mesh used in the simulation using Abaqus pQ. The right figure shows 
contours of concentration obtained using Abaqus for the problem described in subsection 
14.5.11 The white area in the right figure indicates the region in which concentration has a 
negative value. Approximately 37 % of nodes have negative values. The minimum concen- 
tration obtained is —0.0832 (which is -4.16 % off the range of possible values: to 2). Wc 
have taken d\ = 10000, d,2 = 1, and 9 = —ir/6. 
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Figure 2. Numerical /i-convergence: A pictorial description of the boundary value problem. 
The computational domain is a bi-unit square with origin denoted as O. The boundary 
conditions and the volumetric source for the diffusion problem is shown in the left figure. 
The (Dirichlet) boundary conditions and the specific body force for the deformation problem 
is shown in the right figure. The Dirichlet boundary conditions u p (x) are prescribed by 
directly evaluating the expressions for the displacement given in equation (|37p . 
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Figure 3. Numerical /i-convergence: This figure shows the typical three-node triangular 
and four-node quadrilateral meshes used in the numerical convergence analysis. Both these 
meshes have 21 nodes along each side. For the numerical convergence analysis, these meshes 
are subdivided in a hierarchical manner to obtain a scries of meshes. 




Figure 4. Numerical /i-convergence: Comparison of concentration profile from analytical 
solution (left) to that of the numerical study (right) using three-node triangular mesh shown 
in Figure |31 
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Figure 5. Numerical /i-convergence: Convergence of the proposed computational frame- 
work with respect to the concentration field is illustrated in this figure. We show the 
convergence in L 2 -norm and iJ 1 -seminorm for three-node triangular (left) and four-node 
quadrilateral (right) finite elements. 
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Figure 6. Numerical /i-convergence: Convergence of the proposed computational frame- 
work with respect to the displacement field is illustrated in this figure. We show the 
convergence in L 2 -norm and i^-seminorm for three-node triangular (left) and four-node 
quadrilateral (right) finite elements. 
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Figure 7. Cantilever beam with edge shear: A pictorial description of the boundary value 
problem with origin at 'O'. The boundary conditions and the specific body force for the 
deformation problem is shown in the top figure and the boundary conditions and the volu- 
metric source for the diffusion problem is shown in the bottom figure. 
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Figure 8. Cantilever beam with edge shear: Four-node quadrilateral mesh (top figure) 
used in the numerical study and comparison of the concentration profiles (in the order of 
precedence) for three different values of $s being equal to 5, 10 and 20 
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Figure 9. Simply supported beam under self- weight: A pictorial description of the bound- 
ary value problem with origin at 'O'. The boundary conditions and the specific body force 
for the deformation problem is shown in the top figure and the boundary conditions and the 
volumetric source for the diffusion problem is shown in the bottom figure. 
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Figure 10. Simply supported beam under self- weight: Four-node quadrilateral mesh (top 
figure) used in the numerical study and comparison of the concentration profiles (in the 
order of precedence) for three different values of rjs being equal to 1, 10 3 and 2 x 10 4 
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Figure 11. Fixed beam under self- weight: A pictorial description of the boundary value 
problem with origin at 'O'. The boundary conditions and the specific body force for the 
deformation problem is shown in the top figure and the boundary conditions and the volu- 
metric source for the diffusion problem is shown in the bottom figure. 
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Figure 12. Fixed beam under self- weight: Four- node quadrilateral mesh (top figure) used 
in the numerical study and comparison of the concentration profiles (in the order of prece- 
dence) for three different values of $t being equal to 1, 5 and 7. 
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Figure 13. Fixed beam under self- weight: The concentration profiles obtained using the 
diffusivity model (|45|) as outlined in Reference [25] . 
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Figure 14. Plate with a hole under self-weight: A pictorial description of the dimensions 
for a plate with a hole in the left figure and three-node triangular unstructured mesh used in 
our numerical study in right figure. The origin is located at 'O' and the vertex 'H' of the hole 
at (|, |). Displacements at the boundary of the hole and traction at the boundary of the 
plate are zero. In the region between the plate and the hole, the body force is b(x) = — 10 e y 
and the volumetric source is /(x) = 0. The concentration at the boundary of the hole is 1 
and at the boundary of the plate is 0. 
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Figure 15. Plate with a hole under self-weight: This figure compares concentration profiles 
from the standard Galerkin formulation (left figure) and the non-negative formulation (right 
figure). The white area in the left figure indicates the regions in which concentration is 
negative. The violations also include that are of order machine precision. 




Figure 16. Plate with a hole under self-weight: This figure compares degradation profiles 
from the standard Galerkin formulation (left figure) and the non-negative formulation (right 
figure). The red area indicates the regions in which the material is degrading and the green 
area indicates the regions in which the material is healing (which is because unphysical 
negative value for the concentration). In this figure we do not consider negative values of 
order machine precision as violations. 
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Figure 17. Beam with three holes under self- weight: A pictorial description of the dimen- 
sions for a beam with three square holes. The dimensions of each hole are 0.4 x 0.4. The 
origin is located at 'O' and the vertices 'Hi', 'H 2 ' and 'H 3 ' of the holes 'A, B, C are at 
(1.8, 0.3), (4.8, 0.3) and (7.8, 0.3). Displacements at the boundary of the holes and traction 
at the boundary of the beam are zero. In the region between the beam and the holes, the 
body force is b(x) = — 10 and the volumetric source is /(x) = 0. The concentration at 
the boundary of the holes is 1 and at the boundary of the beam is 0. 
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Figure 18. Beam with three holes under self- weight: Comparison of concentration profile 
from standard Galcrkin formulation (middle figure) to that of the non-negative formulation 
(bottom figure). Three- node triangular unstructured mesh used in the numerical study is 
shown in the top figure. The white area in the middle figure indicates the region in which 
concentration has a negative value. In this figure, the white area also includes the regions 
with concentration of order machine precision. 
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Figure 19. Beam with three holes under self- weight: Comparison of degradation profile 
from standard Galerkin formulation (top figure) to that of the non-negative formulation 
(bottom figure). The red area indicates material is degrading and the green area indicates 
the material is healing (because of unphysical negative value for the concentration). In 
this figure we do not consider the negative values that are of order machine precision as 
violations. 
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Figure 20. Convergence of staggered coupling algorithm: We have plotted In[||cW — 
c^ _1 '||] with respect to iteration number for various problems. The stopping criterion in 
the staggered coupling algorithm is ||cW — c( l_1 )|| < e^OL- As one can see f rom the figure, 
the staggered coupling algorithm converges for all the chosen problems. Also, note that the 
algorithm need not converge monotonically, which is evident in the case of plate with a hole 
test problem. 

Maruti Kumar Mudunuru, Graduate Student, Department of Mechanical Engineering, Texas A&M 
University, College Station, Texas 77843. 
E-mail address: maruti.iitm@neo.tamu.edu 

Correspondence to: Dr. Kalyana Babu Nakshatrala, Department of Mechanical Engineering, 216 
Engineering/Physics Building, Texas A&M University, College Station, Texas 77843. TEL: + l-979- 
845-1292 

E-mail address: knakshatrala@tamu.edu 



.39 



10 



G 

o 



I 

> 



Cantilever 



o 

^ ■ Fixed 



▲ 



Simply supp. 



4 5 6 7 8 9 10 11 12 
Staggered iterations 



I 

< 



22 
20 
18 
16 
14 



12 





* -•-Plate with a hole 




g' p ■» Beam with three holes 

J / /•. 7, 




\\ 


S , ' ! * 

* 




»* ( l i 1 

\: n / \ / 






V V i • / 






■ 






■' 









5 10 15 

Staggered iterations 



20 



Figure 21. Convergence of the active-set strategy: This figure shows the number of itera- 
tions taken by the active-set strategy for each iteration in the staggered coupling algorithm 
for various test problems. 
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